﻿#region Copyright information
// 
// Copyright © 2005-2013 Yongkee Cho. All rights reserved.
// 
// This code is a part of the Biological Object Library and governed under the terms of the
// GNU Lesser General  Public License (LGPL) version 2.1 which accompanies this distribution.
// For more information on the LGPL, please visit http://bol.codeplex.com/license.
// 
// - Filename: InverseGaussianDistribution.cs
// - Author: Yongkee Cho
// - Email: yongkeecho@gmail.com
// - Date Created: 2012-09-06 11:39 AM
// - Last Modified: 2013-01-25 3:59 PM
// 
#endregion
using System;
using System.Collections.Generic;
using System.Linq;
using BOL.Maths.Functions;

namespace BOL.Maths.Distributions
{
    /// <summary>
    /// Represents an inverse Gaussian distribution with mean mu and scale parameter lambda]..
    /// <see>Mathematica: ref/NormalDistribution</see>
    /// <see>http://en.wikipedia.org/wiki/Inverse_gaussian_distribution</see>
    /// </summary>
    public sealed class InverseGaussianDistribution : IUnivariateDistribution<double>, IEquatable<InverseGaussianDistribution>
    {
        #region Private variables

        private double _mu; // mean
        private double _lambda; // scale parameter

        #endregion

        #region Public properties

        /// <summary>Gets or sets the mean of an inverse Gaussian distribution.</summary>
        public double Mu
        {
            get { return _mu; }
            set
            {
                if (value <= 0)
                    throw new ArgumentOutOfRangeException("value", "mu is expected to be positive");
                
                _mu = value;
            }
        }

        /// <summary>Gets or sets the scale parameter of an inverse Gaussian distribution.</summary>
        public double Lambda
        {
            get { return _lambda; }
            set
            {
                if (value <= 0)
                    throw new ArgumentOutOfRangeException("value", "lambda is expected to be positive"); 
                
                _lambda = value;
            }
        }

        /// <summary>Gets the range of an inverse Gaussian distribution.</summary>
        public IRange<double> Domain { get { return new Range<double>(0.0, double.PositiveInfinity); } }

        /// <summary>Gets the mean of an inverse Gaussian distribution.</summary>
        public double Mean { get { return _mu; } }

        /// <summary>Gets the median of an inverse Gaussian distribution.</summary>
        public double Median { get { throw new Exception("Median is undefined for inverse Gaussian distribution."); } }

        /// <summary>Gets the mode of an inverse Gaussian distribution.</summary>
        public double Mode { get { return _mu * (Math.Sqrt(1 + 9 * _mu * _mu / (4 * _lambda * _lambda)) - 3 * _mu / (2 * _lambda)); } }

        /// <summary>Gets the variance of an inverse Gaussian distribution.</summary>
        public double Variance { get { return Math.Pow(_mu, 3) / _lambda; } }

        /// <summary>Gets the skewness of an inverse Gaussian distribution.</summary>
        public double Skewness { get { return 3 * Math.Sqrt(_mu / _lambda); } }

        /// <summary>Gets the kurtosis of an inverse Gaussian distribution.</summary>
        public double Kurtosis { get { return 15 * _mu / _lambda; } }

        /// <summary>Gets the entropy of an inverse Gaussian distribution.</summary>
        public double Entropy { get { throw new Exception("Entropy is undefined for inverse Gaussian distribution."); } }

        #endregion

        #region Constructor

        /// <summary>
        /// Instantiates an inverse Gaussian distribution.
        /// </summary>
        /// <param name="mu">mean</param>
        /// <param name="lambda">_lambda</param>
        public InverseGaussianDistribution(double mu, double lambda)
        {
            if (mu <= 0)
                throw new ArgumentOutOfRangeException("mu", "mu is expected to be positive");
            if (lambda <= 0)
                throw new ArgumentOutOfRangeException("lambda", "lambda is expected to be positive");

            _mu = mu;
            _lambda = lambda;
        }

        #endregion

        #region ICloneable implementation

        public InverseGaussianDistribution Clone()
        {
            return new InverseGaussianDistribution(_mu, _lambda);
        }

        object ICloneable.Clone()
        {
            return Clone();
        }
        
        #endregion

        #region Public methods

        /// <summary>
        /// Returns the inverse Gaussian density given value.
        /// </summary>
        /// <param name="value"></param>
        /// <returns>inverse Gaussian density</returns>
        public double Pdf(double value)
        {
            if (value <= 0)
                throw new ArgumentOutOfRangeException("value", "value is expected to be a positive real number.");

            return Math.Exp(-0.5 * _lambda * Math.Pow(value - _mu, 2) / value / Math.Pow(_mu, 2)) * Math.Sqrt(_lambda / Math.Pow(value, 3)) / Math.Sqrt(2 * Math.PI);
        }

        /// <summary>
        /// Returns the cummulative inverse Gaussian density given value.
        /// </summary>
        /// <param name="value"></param>
        /// <returns>cummulative inverse Gaussian density</returns>
        public double Cdf(double value)
        {
            if (value > 0)
            {
                var a = Math.Sqrt(_lambda / value);
                var b = value / _mu;
                var c = Math.Sqrt(2);
                var d = Math.Exp(2 * _lambda / _mu);

                return 0.5 * ((1 + DistributionFunctions.Erf(a * (b - 1) / c)) + d * (1 - DistributionFunctions.Erf(a * (b + 1) / c)));
            }

            return 0;
        }

        /// <summary>
        /// Returns the value of the inverse Gaussian distribution for given probability.
        /// </summary>
        /// <param name="p">probability</param>
        /// <returns></returns>
        public double Quantile(double p)
        {
            throw new NotImplementedException();
        }

        /// <summary>
        /// Returns a random value from the inverse Gaussian distribution.
        /// <see>http://en.wikipedia.org/wiki/Inverse_gaussian_distribution</see>
        /// </summary>
        /// <returns>random value</returns>
        public double Sample(Random r)
        {
            var nu = new GaussianDistribution(0, 1).Sample(r);
            var mu2 = _mu * _mu;
            var y = nu * nu;
            var x = _mu + mu2 * y * .5 / _lambda - _mu * .5 * Math.Sqrt(4 * _mu * _lambda * y + mu2 * y * y) / _lambda;
            var z = r.NextDouble();

            return z < (_mu / (_mu + x)) ? x : mu2 / x;
        }

        /// <summary>
        /// Estimates parameters of the inverse Gaussian distribution using maximum likelihood.
        /// <see>http://en.wikipedia.org/wiki/Inverse_gaussian_distribution#Maximum_likelihood</see>
        /// </summary>
        public void MaximumLikelihoodEstimate(IEnumerable<double> source)
        {
            if (source == null)
                throw new ArgumentNullException("source");

            var list = source.ToList();
            Mu = list.Average();
            Lambda = list.Count / list.Sum(xi => 1 / xi - 1 / Mu);
        }

        /// <summary>
        /// Estimates parameters of the inverse Gaussian distribution using Bayesian methods.
        /// </summary>
        public void BayesianEstimate(IEnumerable<double> source, params dynamic[] priorParameters)
        {
            throw new NotImplementedException();
        }

        #endregion

        #region IEquatable<InverseGaussianDistribution> implementation

        public bool Equals(InverseGaussianDistribution other)
        {
            return _mu.Equals(other._mu) & _lambda.Equals(other._lambda);
        }

        #endregion

        #region Object overriden

        public override int GetHashCode()
        {
            return Mu.GetHashCode() ^ Lambda.GetHashCode();
        }

        public override bool Equals(object other)
        {
            if (other == null)
                throw new ArgumentNullException("other");

            if (!(other is InverseGaussianDistribution))
                throw new InvalidCastException("The 'other' argument is not a InverseGaussianDistribution object.");
            
            return Equals(other as InverseGaussianDistribution);
        }

        public override string ToString()
        {
            return String.Format("InverseNormal(mu = {0:G4}, lambda = {1:G4})", _mu, _lambda);
        }

        #endregion
    }
}
